function spect = spectrum_analyser(spectrum, varargin)
%SPECTRUM_ANALYSER Summary of this function goes here
%   Detailed explanation goes here
    p=inputParser;
    isnonneg = @(x) x > 0;
    spectrumTest = @(x) all(ismember({'frequency', 'value'}, fields(x)));
    freqTest = @(x) x >= min(spectrum.frequency) && x <= max(spectrum.frequency);
    p.addRequired('spectrum', spectrumTest);
    p.addParameter('nbins', 100, isnonneg);
    p.addParameter('fmin', min(spectrum.frequency), freqTest);
    p.addParameter('fmax', max(spectrum.frequency), freqTest);
    p.addParameter('isPlot', true, @islogical);
    p.addParameter('NPeaks', 1, isnonneg);
    p.addParameter('MinPeakHeight',-Inf, @isnumeric);
    p.addParameter('MinPeakProminence',1e-10, @isnumeric);
    p.addParameter('NENBW', 1.5, isnonneg);
    p.addParameter('figure_hdl', 'new');
    p.parse(spectrum, varargin{:});
        
    freq_all = p.Results.spectrum.frequency; 
    
    df = freq_all(2) - freq_all(1);
    spect.f_resolution = df;
    spect.nenbw = p.Results.NENBW;
    spect.enbw = p.Results.NENBW*df;

    val_all = p.Results.spectrum.value/sqrt(spect.enbw);

    spect.all.frequency = freq_all;
    spect.all.value = val_all;
    
    %%
    if strcmp(p.Results.figure_hdl, 'new') && p.Results.isPlot
        figure; hdl = axes;
    else
        hdl = p.Results.figure_hdl;
    end
    %% noise backgournd
    
    
    freq_sel_idx = intersect(find(freq_all < p.Results.fmax), ...
                             find(freq_all > p.Results.fmin) );
                         
    freq_sel = freq_all(freq_sel_idx);
    val_sel = val_all(freq_sel_idx);
    
    nbins=p.Results.nbins;
    

    [count, edge] = histcounts(log10(val_sel), nbins);
    [~,idx] = max(count);
    l=edge(idx);r=edge(idx+1);

    bg_idx = intersect(find(log10(val_sel)< r), find(log10(val_sel)>l));
    
    
    spect.background.frequency = freq_sel(bg_idx);
    spect.background.value = val_sel(bg_idx);
    spect.background.mean_level = mean(spect.background.value);
    spect.background.level_std = std(spect.background.value);
    
    if p.Results.MinPeakHeight == -Inf
        minPkH = 5.0*spect.background.level_std + spect.background.mean_level;
    else
        minPkH = p.Results.MinPeakHeight;
    end
    minPkPro = p.Results.MinPeakProminence;

    %% peak
    [pks,locs,w,pr] = findpeaks(log10(val_sel), ...
                                'NPeaks', p.Results.NPeaks, ...
                                'MinPeakHeight', log10(minPkH), ...
                                'MinPeakProminence', minPkPro );
    spect.peaks.value = 10.^pks;
    spect.peaks.loc = locs;
    spect.peaks.width = w;
    spect.peaks.prominence = pr;
    
    spect.peaks.data=cell(1, p.Results.NPeaks);
    for k=1:p.Results.NPeaks
        z = 5*ceil(w(k)/df);
        l=max(1, locs(k)-z); r=min(locs(k)+z, length(freq_sel) ); 
        peak_k.frequency = freq_sel(l:r);
        peak_k.value = val_sel(l:r);
        spect.peaks.data{k} = peak_k;
    end

    if p.Results.isPlot
        semilogy(hdl, ...%freq_all, val_all, '--r', ...
                 freq_sel, val_sel, '.-r', ...
                 freq_sel(locs), val_sel(locs), 'dk', ...
                 ...%spect.background.frequency, spect.background.value, '.r', ...
                 minmax(freq_all), [spect.background.mean_level, spect.background.mean_level], 'k--')
        xlim(hdl, minmax(freq_sel)); grid on;
        text(hdl, double(min(freq_sel)), ...
                  double(spect.background.mean_level + 10*spect.background.level_std), ...
                  sprintf('Noise floor\n%3.2e V\\cdot Hz^{-1/2}', spect.background.mean_level), 'FontSize', 15);
        for k=1:length(spect.peaks.data)
            text(hdl, double(freq_sel(locs(k))+3*df), ...
                      double(val_sel(locs(k))), ...
                      sprintf('%3.2e V_{pk}',val_sel(locs(k))*sqrt(spect.enbw) ), 'FontSize', 15);
        end
        title(sprintf('f_{resolution}=%3.2f Hz, ENBW = %3.2f Hz, NENBW=%2.1f', spect.f_resolution, spect.enbw, spect.nenbw));
            
    end
end

